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1.0 INTRODUCTION 


The objective of the study described in this report is to investi- 
gate the influence of physical constraints which may be approximately satis- 
fied by the Earth's liquid core on models of the geomagnetic main field and 
its secular variation. A previous report (Estes, 1984) describes the metho- 
dology used to incorporate non-linear equations of constraint into the main 
field model, and will be referred to as Report I. This report describes the 
application of that methodology to the GSFC 12/83 field model (Langel and 
Estes, 1985a) to test the frozen-flux hypothesis and the usefulness of incor- 
porating magnetohydrodynamic constraints for obtaining improved geomagnetic 
field models. 

The geomagnetic field is a highly complex phenomenon with signifi- 
cant spatial and temporal variations. Global satellite surveys such as P0G0 
and MAGSAT have led to important advances in our understanding of this global 
geophysical phenomenon and its applications in geology, tectonophysics, navi- 
gation, communications, and many other fields. The geomagnetic field is com- 
posed of three major components: 

• The core field is the major part of the geomagnetic field 
and is believed to be caused by fluid motion of conducting 
material in the Earth's outer core. The field variations 
have large amplitudes, broad spatial features (i.e., wave- 
lengths on the order of 1000 km and longer at the earth's 
surface), and significant low-frequency temporal 
variation. 

• The crustal anomaly field forms a relatively small part of 
the geomagnetic field, but is significant since it is 
caused by geologic structure and the state of crustal 
minerals. The spatial variations of the crustal field 
occur with characteristic wavelengths ranging up to a few 
thousand km. 
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• The external field is caused by electrical currents in the 
ionosphere and magnetosphere. Magnetic field variations 
due to the external field occur on all spatial and 
temporal scales. The external field accounts for the most 
complex and dynamic part of the geomagnetic field. 

Models of the main, or core field are required for a number of uses, 
including the production of declination charts for navigation, the removal of 
background fields from aeromagnetic and ship-borne magnetic surveys, satellite 
attitude control and determination, the probing of various fluid motions in 
the Earth's core, and calculations of field lines and conjugate points for 
ionospheric and magnetospheric studies. The geomagnetic secular variation has 
been demonstrated to be correlated with decade fluctuations in the Earth's 
rotation, and improved models could lead to better understanding of the nature 
of the core/mantle coupling. Moreover, knowledge of the secular variation of 
the core field is crucial in tying together regional aeromagnetic and marine 
surveys taken at different epochs. The secular variation is one of the few 
sources of information regarding the Earth's core, where the internal field 
originates. 

As indicated above, the anomaly portion of the field arises from 
variably magnetized rocks lying between the Earth's surface and the Curie 
isotherm (the "magnetic crust"), and is the part of the magnetic field having 
direct geologic significance. The magnetic anomaly fields, which are essen- 
tially time independent for our purpose, have long been used to study near- 
surface geology and structure. They have had an important role in the initial 
evidence for sea-floor spreading that led to our present unifying concepts of 
plate tectonics and they are routinely used in resource exploration. Models 
of the sources of magnetic anomalies are the basis for geologic interpreta- 
tion. 


The modeling of the main and anomalous fields with high resolution 
is of major importance for earth science applications. In order to achieve 
this resolution, techniques for developing core and crustal models must be 
improved. Separating the core field from the crustal field is a major 
problem. Conventional approaches to representing the main field by a 
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truncated spherical harmonic series suffer from aliasing due to the crustal 
field above approximately degree and order twelve (Langel and Estes, 1982). 
This problem is described in detail by Langel (1985). A consequence of this 
is that extrapolation of main field models to the core/mantle boundary (CMB) 
is severly contaminated by the attenuated higher harmonic coefficients which 
are aliased by crustal influences. New approaches to this problem by Gubbins 
(1983, 1984) and Gubbins and Bloxham (1985) using a priori spectral assump- 
tions and by Shure et. al. (1982) using harmonic splines minimizing a smooth- 
ness norm for the field at the CMB are important contributions. Currently, 
Langel and Estes (1985b) are investigating the separability problem through a 
correlated measurement noise weight matrix approach. 

In addition to the ambiguity inherent in the separation process 
which imposes a fundamental limitation on accuracy, there is further difficul- 
ty in resolving the non-static nature of the main field. The main field 
changes in time on scales ranging from seconds to hundreds of thousands of 
years. Magnetic polarity reversals, for example, can take place over a period 
of a few thousand years, while tens to hundreds of thousands of years elapse 
between reversals. The interest of our study, however, is on magnetic varia- 
tions over the scale of several decades and less. The secular variation (SV), 
even over those time scales, is less well known than the main field. This is 
due primarily to data availability. While the dense global vector magnetic 
survey provided by Magsat (Langel et. al, 1980) has permitted an excellent 
global "snapshot view" of the main field at epoch 1980, secular variation 
models depend principally on the data from the global network of permanent 
magnetic observatories which are few in number and poorly distributed. An 
additional difficulty is that magnetic features of smaller spatial scale are 
probably capable of more rapid temporal variation so that higher degree terms 
of the spherical harmonic series representation of the secular variation may 
be more important and convergence of the series weaker than that for the 
constant part of the main field. The result is that with current methods of 
modeling the main field and its secular variation even the large scale compon- 
ents of the geomagnetic field cannot be predicted into the future by more than 
a few years beyond the data span which defined the model. Two approaches are 
possible to improve the situation; to incorporate some statistical information 
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describing the SV process into the model (Gibbs and Estes, 1982, 1984), or to 
insert physical constraints provided by the electromagnetism and fluid dynam- 
ics of the source region, earth's outer core. 

The work presented in this report describes the research accomplish- 
ed in the second year of a multi-year study of the influence of magnetohydro- 
dynamic constraints from the core on geomagnetic models, which is a joint 
effort between Edward Benton of the University of Colorado and Business and 
Technological Systems, Inc. Some theoretical aspects of the problem are pre- 
sented by Benton (1984, 1985) in separate reports. As described in Report I, 
the BTS contribution in the first year was devoted to establishing appropriate 
geomagnetic data bases, establishing algorithms for incorporating non-linear 
constraints into the estimation procedure, and developing and verifying 
software based on these algorithms for computation. The demonstration that 
geomagnetic observations are consistent with the field being frozen to the 
core fluid has great importance for computing properties of the fluid core 
velocities, and for constructing improved predictive models. Gubbins (1984) 
and Bloxham and Gubbins (1985) have investigated the consistency of the 
observed data with frozen-flux constraints by two methods; they developed 
models that were forced to satisfy the flux constraints and then examined for 
acceptability the resulting misfit to the data used in the model, and they 
computed the flux integrals from unconstrained models at different epochs and 
compared their differences with formal error estimates from the models at the 
core surface. Our approach in this study uses the methodology of Report I to 
produce a constrained geomagnetic model based on GSFC (12/83) using the 
constraint equations as additional "data" with an associated uncertainty, and 
compares the misfit errors for data beyond the time interval defining the 
model (prediction errors) with those of the unconstrained model. 


4 


Bvsiness .ind Technological Systems, Inc 


2.0 PHYSICAL CONSTRAINTS FOR MAIN FIELD MODELS 


In the earth's core the geomagnetic field satisfies the induction 

equation 


B_ = V * {V*b.) + nV 2 JJ (1) 

where is the vector field, is the flow velocity and n is the 
magnetic diffusivity. On the time scale of a few decades, it has been shown 
to be a good approximation to ignore magnetic diffusion in the core (Backus, 
1968), so that the core fluid is assumed to behave like a perfect conductor 
with the magnetic flux frozen into the fluid motion. Moreover, it is a 
reasonable first order approximation to treat the mantle as a perfect 
insulator. Horizontal components of J3 may be discontinuous across the CMB, 
but the radial component, B r , must be continuous. Bondi and Gold (1950) 
introduced the scalar relation for the absolute unsigned magnetic flux, or 
"pole-strength" 


M(r, t) = Sf /J |B r (r, 9, A, t) | r 2 sine d9dx (2) 

as a measure of the number of field lines which cross the spherical surface of 
radius r at time t . They noted that the number of field lines crossing 
the top surface of a perfectly conducting core (CMB) is a constant of the 
motion. Then 


M(b, t) = M q = constant (3) 

is a scalar constraint for the frozen-flux hypothesis, where b is the radius 
of the CMB. Voorhies and Benton (1982) show that this constraint appears to 
be well satisfied, which provides support for the perfect conducting 
core/insulating mantle hypothesis over short time scales. Moreover, Backus 
(1968) has shown that the total absolute flux linking the CMB is the sum of 
contributions arising from separate "null-flux" patches, which are regions on 
the CM8 bounded by contours of zero radial field. Denoting the contour of the 
i^ 1 null -flux curve as r. , then 
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M . (b, t) = / L B (b, 0, x, t) b 2 sine dedx (4) 

1 l r 

= constant 

provide additional scalar integral constraints. Unfortunately, the null 

curves r. are non-linear functions of the geomagnetic Gauss coefficients 

and, moreover, depend strongly on the truncation level of the spherical 

harmonic series for . 

r 

Benton (1984, 1985) has pursued additional physical constraints on 
secular variation arising from the coupling of non-dissipative fluid dynamics 
to the frozen-flux electromagnetism of earth's core. He shows that the signed 
flux through any patch which is bounded by a curve which moves with the fluid 
should be conserved in the limit of perfect core conductivity, and examines 
possible candidates. He demonstrates that under some assumptions, for 
example, lines of absolute vorticity move with the fluid, and the contours on 
which the vertical absolute vorticity vanishes move with the fluid. A 
consequence of this is that under the stated assumptions (spherical, perfectly 
conducting in viscid core free of vertical Lorentz torque at the CMB) the 
geomagnetic flux through the northern geographic hemisphere on the core-mantle 
boundary (CMB) should be constant. With the assumption that the mantle is an 
insulator, the resulting constraint equation 

L E = fo" r 0 /2 B r (b, 0, *, t) b 2 sin© d0d<|. 

N n+2 n . 

= n 9 n (t) P H ( 7 ) = COnSt * {5) 

is linear in the Gauss coefficients g nm . Here a is the Earth's radius, b 
is the radius of the CMB, and P* is the Schmidt normalized Legendre poly- 
nomial of degree n and order 1. Denoting 
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L e = l F n 9p(t) = constant 
m=l 


we have 


c _ _ ( n ~2 ) c 

F n " HvTT ( b J F n-2 


where Fg = 0 and = 1. Note that this represents a constraint equation 
for secular variation 


N n 

L E = l F n 5„(t) = 0 (6) 

L n=l " n 

which may be readily tested against existing models or easily incorporated as 
a linear constraint into a model. 

Linder the assumptions which established the flux linking the 
northern geographic hemisphere as a constant, additional flux constraints may 
be established for areas bounded by curves defined by the intersection of the 
geographic equator with the null curve along the geomagnetic equator on the 
CMB. These constraints, like those due to the frozen flux through the null 
curves, are non-linear functions of the Gauss coefficients and must be 
implemented numerically. 
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3.0 BAYESIAN LEAST SQUARES PARAMETER ESTIMATION 


The Bayesian parameter estimation algorithm provides a methodology 
for including a priori statistical information on the parameter space in 
obtaining least squares solutions (Luenberger, 1973). Let x denote the 
parameter vector, y the measurement vector and v represent the random 
noise vector with zero mean and covariance matrix R . Then the observation 
equation is 


y = F ( x ) + v . (7) 

In the linear case, assume F ( x) = Ax so that 

y = Ax + v . (8) 

Then if x g denotes an a priori estimate of the parameters with a priori 
covariance matrix , the estimate for the parameter state vector x is 

a 

x = (A T R _1 A + fi a ’ 1 ) 1 [A^'V + o’ 1 x a ] . (9) 

This result may be obtained by minimizing the norm 

J,c(x) = (y - Ax) T R _1 (y - Ax) + (x - x a ) T fl^fx - x) (10) 

lo ad a 

with respect to x , so that it is clear that the a priori information is 
included in the formalism as additional data. 

For the non-linear problem, an iterative approach is required. 
Linearizing about a nominal solution Xq , we have 

F(x) - F(Xq) = A(xq)(x - x Q ) + .... . (11) 
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s t 

A Gauss iteration procedure yields the equation for the (n+1) approximation 
to the solution estimate where 


x n+l = x n + (5x n+l 


( 12 ) 


6x n +i = (A T (x n ) R _1 A(x n ) 



[A (x n ) R 6y n + ^a^[x a - 



and 

Sy„ = y - A(x ) x n 
n n n 

The rate of convergence of the procedure depends on how good the nominal guess 
x Q is, and on the non-linear character of the function F(x) . For highly 
non-linear problems, other techniques utilizing higher order derivatives may 
be required. 

The implementation of constraints may be accomplished within the 
estimator formalism by Lagrange multiplier techniques (see Gubbins, 1984) or 
by considering the constraint equations to be data. If a non-linear 
constraint equation is represented as 

G(x) = g (13) 

then the (n+l) st iteration relation is 

8 Vl ■ [ A T ( X n ) R _1 A(x n ) + a " 1 * C T (x n ) R' 1 C(x n )]" 1 (14) 

• ( A T (; n > R _1 «y„ ♦ - x n ) * C T (x n ) R^gJ 


where 


C(x ) = — 
9 x 


( 15 ) 


<sg n - 


- C(x n )x n 
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The weight matrix R * reflects the stiffness of the constraint, where R 

c c 

is the “covariance" matrix of the "observed" constraint. This error measure 
may be chosen to represent the estimated numerical precision lost in the 
computing process, or to represent modeling errors inherent in the constraint. 
This approach has also been taken by Bloxham and Gubbins (1985). 

In the absence of constraints, the least squares estimation process 
converges rapidly (in usually no more than two iterations) when the time 
dependence of Gauss coefficients is a power series. In that representation, 
the observations of magnetic components X, Y, Z , are linear functions of the 
Gauss coefficients whereas the scalar field intensity, horizontal field, 
declination, and inclination are nonlinear functions of those parameters. 

While the measurements will be nonlinear functions of the parameters, they are 
at least, still represented in analytic form. In contrast, some of the 
constraints are not available in analytic form. For example, the null-curves 
r . on which the vertical magnetic field vanishes are extremely complex 
functions of the model parameters and must be determined numerically. As a 
result numerical methods must be used to compute the matrix 

C T (x )C ( x ) 

-n -n 


for each iteration, n = 0, 1, 2, 


until convergence is attained. 
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4.0 EVALUATION OF A FLUX CONSTRAINED GEOMAGNETIC FIELD MODEL 


To investigate the influence of physical constraints which may be 
satisfied by Earth's liquid core on geomagnetic field models, we apply the 
non-linear algorithms described in Section 3.0 to a variation of the GSFC 
12/83 model (Langel and Estes, 1985a). This model utilized quiet Magsat data 
and observatory annual means for a selected set of observatories for the years 
1977 through 1982. The internal Gauss coefficients were expanded to degree 13 
for the main terms and degree 10 for the secular variation terms. In addition 
to the internal Gauss coefficients, GSFC 12/83 solved for an external field of 
spherical harmonic degree one, and vector biases for each observatory. More- 
over, parameters for linear D t variation in g^ Q and the external coeffi- 
cients are included in the solution. Internal coefficients which were not 
determined in the inversion of the data to the 95 % confidence level were 
forced to be zero in the final solution. Because the incorporation of dynamic 
core constraints could improve the observability of such suppressed coeffi- 
cients, the GSFC 12/83 model was recomputed with all coefficients retained in 
the solution and all partial derivatives included in the normal matrix 
A T R -1 A. This model is designated GSFC 9/84-0 and is summarized in Table 1. 

The GSFC 9/84-0 model, with the utilization of high quality Magsat 
data, very accurately represents the geomagnetic field at and above the 
Earth's surface at epoch 1980.0. Figure 1 shows the radial component of the 
field from GSFC 9/84-0 at the core-mantle boundary (CMB) at 1980.0, using all 
Gauss coefficients to degree and order thirteen. There are eleven separate 
null flux curves resulting from this model evaluation. Because the higher 
order Gauss terms, which are less well determined, become more dominant in the 
field evaluation at the CMB, the resulting computations are subject to large 
uncertainty. Benton, et.al. (1982) and Report I examined the effect of 
truncation level on geomagnetic properties at the CMB, and demonstrated 
considerable variation of null-curve location and null-patch flux as a 
function of degree and order. For this study, we have selected the five 
null-curves indicated in Figure 1 for which to impose flux constraints because 
these curves qualitatively demonstrated stability of form and location from 
truncation level nine through thirteen. Moreover, the consideration of only 
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five null -patches meets the objectives of the study while keeping the 
computational burden within reasonable bounds. 


The methodology for computing the null -patch flux and the partial 
derivatives of the null -patch flux with respect to the model parameters (Gauss 
coefficients) in equation (15) is described in Report I. A numerical 
derivative is taken for each parameter using a 1° x 1° grid overlay of the 
null-patch and varying the nominal Gauss coefficient value. In this 
approximation the field is computed at the center of each 1° x 1° cell and 
assumed to be constant within the cell. Note that the flux depends on the 
constant and secular variation Gauss coefficients both through the field 
evaluation and the null-curve boundary, r.. 

For each of the five null -curves selected, we impose the "data" 

(b, t) = constant 

from equation (4) at the times 1977.5, 1980.0, and 1982.5. Because the GSFC 

9/84-0 model accurately represents the field at epoch 1980.0, we used the 

1980.0 computed flux value for each null-patch to determine the constants. We 

are thus adding fifteen non-linear observations to the GSFC 9/84-0 data set 

which will influence the solution in such a way as to attempt to maintain the 

constancy of the flux through the five null -curves at the given epochs. The 

elements of the diagonal weight matrix R * were set to represent an error 
-5 ^ 

measure of 10 MWb for the "observed" flux. This value was chosen 
primarily by numerical experimentation to provide a stable inversion of the 
normal matrix. 

The solution was taken through three iterations using equation (14) 
starting with the GSFC 9/84-0 model as the nominal starting vector for 
iteration one. No a priori information was assumed in the solution for the 
nominal model, so that n * is zero in equation (14). Each iteration was 

a 

performed in two stages. The first stage processed the Magsat and observatory 
data set as published for GSFC 12/83 and accumulated the normal matrix 
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A^(x ) R * A(x ) and the vector A T (x ) R * 6y . From this information an 
n n n n 

"unconstrained" model was generated for the iteration. These models are 
denoted 6SFC 9/84-1, GSFC 9/84-2, and GSFC 9/84-3 for iterations one, two and 
three, respectively. Note that for iteration one, GSFC 9/84-1 is identical to 
GSFC 9/84-0 since the nominal model has already been converged using the GSFC 
12/83 data set. The second stage for each iteration computed the flux and 
flux partial derivatives for the fifteen "observations", formed the matrix 
C T (x n ) R (; 1 C(x n ) and the vector C^"(x n ) R c * <5g n , and added these to the 
first stage quantities to obtain the "constrained" solution for the 
iteration. The constrained models are denoted GSFC 9/84— 1C, GSFC 9/84-2C and 
GSFC 9/84-3C for iterations one, two, and three, respectively. 

The flux at 1977.5, 1980.0 and 1982.5 and flux rate at 1980.0 
through the five null -patches for all of the models through three iterations 
are presented in Table 2. The maximum variation of the flux over the 1977.5 - 
1982.5 interval for each patch is also given, as well as the RSS of the flux 
rate over all five patches. For a perfectly constrained solution, the flux 
rate through each null -patch should be zero. There is a clear trend showing a 
reduced total flux variation for the five year data interval as the solution 
converges. The poorest improvement is with null-curve 2, which has a maximum 
variation of 30 MWb for GSFC 9/84-0 and 19 MWb for GSFC 9/84-2C. Marked 
improvement is shown for null-curves 3, 4 and 5. The RSS of the flux rate at 
1980.0 through all five patches displays a reduction from 16.1 MWb/yr to 5.9 
MWb/yr. 


The cost function being minimized by the estimation procedure is 

J. = 6y T R" 1 6y + 5g T R -1 Sg 
LS J n J n s n c s n 

The first term represents the weighted sum of squares of the data misfit from 
the Magsat and observatory data, which we denote by Q. Table 3 presents the 
misfit to the Magsat and observatory data for the iterated solutions. As 
expected, the misfit increases with the inclusion of constraint "data". While 
the solution has apparently converged with respect to the Magsat data, the 
situation is not as clear with respect to the observatory data which is more 
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strongly coupled with the constraint “data" through the secular variation. 
Table 4 displays the compliance of the iterated solutions to the linear 
vorticity constraint given in equation (6). 

A stern test for the validity of the null -flux approximation and the 
utility of including physical constraints into geomagnetic field models is the 
demonstration of improved model predictive capability beyond the data interval 
defining the model. Tables 5-8 display the standard deviation of fit to 
observatory data for the 91 observatories used in the model on a yearly basis 
for B, X, Y, and Z for data from 1975.5 through 1984.5. The observatory 
biases recovered with the solutions were included in the model evaluations, 
and the standard deviations were computed unweighted. The constrained models, 
particularly for X and Z, show a clearly improved predictive capability 
beyond the interval 1977.5 - 1982.5 in both directions. Within the resolution 
of the plots, the symbols for GSFC 9/84-3 and GSFC 9/84-3C fall directly on 
GSFC 9/84-2 and GSFC 9/84-2C, respectively, and consequently are not 
displayed. The solid vertical line at 1982.5 is to emphasize that only eight 
observatories were used in the solution for that year, so that the statistics 
for 1982.5 primarily represent predictive errors. The statistics for 1984.5 
are based on a sample of only nine observatories and must be considered 
unreliable. The improved model prediction capability demonstrated 
statistically in Tables 5 through 8 has been achieved without changing the 
mathematical model of the secular variation. Only additional "data" have been 
included. This makes a strong statement for the potential usefulness of 
including physical constraints in field modeling. 

The Gauss coefficient differences between GSFC 9/84-0 and GSFC 
9/84-2C are presented in Table 9. As expected, the effect of adding flux 
constraints is manifest almost totally in the secular variation coefficients, 
with the largest influence principally in the sectorial terms. Differences 
between all GSFC 9/84-3C and GSFC 9/84-2C Gauss coefficients are small and 
show no obvious trends. 
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5.0 FUTURE RESEARCH PLANS 

The comparison of the predictive capability of constrained and 
unconstrained models described in Section 4.0 demonstrates a positive improve- 
ment achieved by incorporating physical null flux constraints. This tends to 
support the consistency of the hypothesis that the field is frozen to the core 
fluid with geomagnetic observations. These results, however, are based on 
utilizing a limited number of null flux constraints and a limited time inter- 
val overwhich the Gauss coefficients were modeled as linear in time. To com- 
plete the analysis begun using the GSFC 12/83 model, a converged constrained 
model will be developed which includes the available observatory data to date 
(1977.5 - 1984.5) and a complete set of null flux constraints. This model 
will be compared with IGRF candidate models for 1985.0. 

To more fully examine the influence of physical constraints in 
modeling secular variation, a goal for future study will be to model the main 
field over a longer time interval utilizing additional null flux constraint 
equations, and to examine the predictive error of constrained and unconstrain- 
ed models. Moreover, examination of Benton's (1985) constraints on secular 
variation based on electromagnetism and fluid dynamics considerations offers a 
linear and a number of additional non-linear constraints which should be 
tested for compliance in unconstrained models and implemented into constrained 
models. Alternate representations to polynomials in time for secular varia- 
tion should also be investigated. 
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